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Abstract 

Local grid refinement aims to optimise the relationship between accuracy of the results and number of 
grid nodes. In the context of the finite volume method no single local refinement criterion has been 
globally established as optimum for the selection of the control volumes to subdivide, since it is not easy 
to associate the discretisation error with an easily computable quantity in each control volume. Often the 
grid refinement criterion is based on an estimate of the truncation error in each control volume, because 
the truncation error is a natural measure of the discrepancy between the algebraic finite-volume equations 
and the original differential equations. However, it is not a straightforward task to associate the truncation 
error with the optimum grid density because of the complexity of the relationship between truncation 
and discretisation errors. In the present work several criteria based on a truncation error estimate are 
tested and compared on a regularised lid-driven cavity case at various Reynolds numbers. It is shown 
that criteria where the truncation error is weighted by the volume of the grid cells perform better than 
using just the truncation error as the criterion. Also it is observed that the efficiency of local refinement 
increases with the Reynolds number. The truncation error is estimated by restricting the solution to a 
coarser grid and applying the coarse grid discrete operator. The complication that high truncation error 
develops at grid level interfaces is also investigated and several treatments are tested. 
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1. Introduction 

The finite volume method is a popular method in Computational Fluid Dynamics (CFD): The domain 
is divided into a number of Control Volumes (CVs) using a grid, the differential equations are integrated 
over each CV, and the integrals are approximated by algebraic expressions involving the values of the 
unknown variables at specific discrete locations (e.g. the CV centres). This results in a system of algebraic 
equations whose solution gives approximate values for the unknowns at these locations. In general, the 
finer the grid the smaller the discrepancy between the algebraic and the differential equations, and the 
better the results. On the other hand, using a higher grid density also increases the computational effort 


* Corresponding author. Tel.: +30 24610 56711; Fax.: +30 24610 21730 
Email address: asirakos@uowm.gr, alexandros.syrakos@gmail.com (Alexandros Syrakos) 


Preprint submitted to Journal of Computational Physics 


August 11, 2015 





required. Therefore, for given computational resources, it is desirable to have a way of determining the 
optimal grid density distribution in space which results in the maximum possible accuracy. 

Some methods move the existing nodes of the grid to new locations in order to improve the accuracy 
of the results without altering the number of nodes of the grid (an interesting example is [1], based on 
truncation error). However, the most popular methods are local refinement methods which are more easily 
applicable to more general grids and complex geometries and which locally increase the grid density by 
subdividing selected CVs into smaller CVs. Various grid refinement criteria have been proposed which 
try to identify the CVs whose subdivision would be most beneficial in terms of increase of accuracy, but 
there is no global agreement among researchers as to which criterion is optimal. 

In the context of the finite volume method, refining CVs where the gradient of a selected flow variable 
is large appears to have been popular in the past (e.g. [2, 3]). This simple choice is useful for capturing 
flow features such as fronts or shock waves, but otherwise it is not appropriate; second- or higher-order 
discretisation schemes should be able to capture these gradients no matter how high without special grid 
refinement, as long as the higher-order derivatives of the flow variables are small. Also, such a method 
could potentially not converge to the exact solution of the partial differential equation as convergence 
may require also the refinement of CVs where the gradients are small but the higher-order derivatives are 
large. 

Another popular class includes local refinement criteria which are based on the truncation error (e.g. 
[4-12]). This has a more solid theoretical background since the truncation error acts as the source of the 
discretisation error (see Section 2). An interesting alternative to the use of the truncation error is the 
use of the residual (e.g. [13, 8, 14-19]): While the truncation error is obtained by applying the discrete 
finite volume operator to the exact solution, the residual is obtained by applying the exact differential 
operator to the approximate finite volume solution. These two quantities have similar properties. Yet 
another method is to base the local refinement criterion on an estimate or indicator of the discretisation 
error itself (e.g. [20, 21]), but as discussed in ([16, 18]) this is not efficient since the discretisation error is 
convected and diffused in the domain and may be high in regions which do not really need a higher grid 
density. Grid refinement indirectly reduces the discretisation error by reducing the source of this error 
(which can be viewed to be either the truncation error or the residual) - and so it is more appropriate to 
base the refinement criterion on the source. 

The present study attempts a detailed evaluation of the efficiency of truncation error-based local 
refinement for the Navier-Stokes equations, by testing and comparing various truncation error-based 
criteria on a lid-driven cavity problem. Although the aforementioned studies demonstrate that truncation 
error-guided local refinement can be beneficial, such a detailed study is missing, to the authors’ knowledge. 
A possible explanation for this is that it is difficult to quantify the efficiency (the ratio of the discretisation 
error reduction to the number of control volumes added to the grid), as the Navier - Stokes equations do 
not have analytic solutions except for very simple cases, and so it is difficult to calculate the discretisation 
error accurately. In the present work a variant of the lid-driven cavity test case is selected (by far the 
most popular test case for incompressible flows), and the problem is first solved on very fine structured 
grids for various Reynolds numbers. Using Richardson extrapolation “exact” solutions are obtained, and 
the various local refinement techniques are evaluated by comparing the discretisation error against the 
number of CVs of the locally refined grids in each case. Sections 2 and 3 present some theory concerning 
the truncation error, its use as a refinement criterion, and a method to estimate it a posteriori Section 
4 describes the specific finite volume method used to solve the incompressible Navier - Stokes equations. 
Finally, Sections 6 and 7 present the lid-driven cavity test case and the results. 

2. The truncation error as a grid refinement criterion 

The truncation error quantifies the discrepancy between the differential equation integrated over a CV 
and the approximate algebraic equation which is obtained for that CV via the finite volume method. This 
makes it appear to be the right quantity to use as a refinement criterion, and indeed many researchers 
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have used it as such (e.g. [4-12]), but there are some complexities. At this point it is appropriate to 
include a brief description of the truncation error and of some of its properties that are of interest. 

Suppose a finite volume method is used to solve the (possibly non-linear) partial differential equation 
(PDE) N(u) = b where N is the differential operator, u is the unknown function, and b is a known 
function independent of u. This equation is integrated over each CV and the result is approximated by 
an algebraic equation. For the whole grid one obtains a system of algebraic equations: 

K 

= N h (u h ) + T h ( 1 ) 

p =1 

Both sides of (1) are vectors with K components, where K is the number of CVs in the grid. On the left- 
hand side, each component P is the integral of the image N(u) over CV P divided by the CV volume 5Vtp. 
On the right-hand side the term Nh(uh) is a vector which is obtained by applying the algebraic operator 
Nh to the vector uh of the values of the exact solution function u at the CV centres. The symbol h is 
sometimes used to denote a particular grid, and sometimes used to denote the spacing between successive 
grid lines of that grid - this will be clear from the context. N\ is constructed by discretising the continuous 
differential operator N according to the finite volume methodology using selected discretisation schemes. 
No matter how accurate these discretisation schemes are, it is in general not possible for the term Nh(uh) 
to equal the left hand side, but they will differ by the truncation error r^. The truncation error consists 
of all the terms that were truncated from the Taylor series used in converting the differential equations 
to algebraic form (see e.g. [22]). 

The PDE N(u) = b implies that both sides of (1) should equal the vector b h of the integrals of b 
over each CV of the grid h divided by the corresponding volume. As the truncation error is unknown, 
the finite volume method drops it under the assumption that it is small enough such that its omission 
does not have a significant impact on the results. Therefore, the system Nh(uh) = bh is solved instead of 
Nh(uh) + Th = frfr, where uh is an approximation to the exact solution and their difference is called 
the discretisation error — Uh- 

The truncation error acts as the source of the discretisation error ([22]), which is distributed in space 
according to the same laws described by the differential equation - e.g., convection and diffusion. This 
can be seen from the following equation which is straightforward to derive: 

A r h{uh + e h ) - N h (u h ) = -r h (2) 

In the linear case this simplifies to N^h — ~ T h where it is evident that the discretisation error obeys 
the same PDE with the truncation error acting as the source. In the non-linear case for small enough 
equation (2) becomes N' h 6h & —where N' h is the Jacobian matrix of the operator N^. This makes it a 
reasonable choice to base the refinement criterion on the truncation error. Local grid refinement normally 
leads to a local truncation error reduction, and according to the above discussion this should also reduce 
the discretisation error, provided that Nh is continuous and is small enough. 

On the other hand there are some complexities. The first is that Th is unknown; however, there are 
ways to estimate it, so this is not a major problem in many cases. Also, the relationship between 
and Th is not simple, as can be seen from (2). If the truncation error in any two CVs is the same, then 
this does not necessarily mean that subdivision of either of them will bring about the same reduction in 
discretisation error. Therefore to use the truncation error as a local refinement criterion, its value in each 
CV should be multiplied by a weighting factor which should come from equation (2). One of the simplest 
options is to divide the truncation error in each CV by the corresponding main diagonal coefficient of 
the Nh operator after it is linearised. This is equivalent to estimating by performing a single Jacobi 
iteration on (2) (after linearising Nh) with zero initial guess, and then using this estimate of as the 

refinement criterion. This sort of weighting has been used e.g. by [13, 14]. 

There exists another difficulty in using the truncation error as the refinement criterion. When the 
grid is smooth then the rate of convergence of the truncation and discretisation errors is the same - 


- 2 — / N(u)dQ 
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e.g. 0(h 2 ) for second order methods. But when the grid is not smooth, e.g. it exhibits significant 
skewness or other geometric irregularities, then it is usual that the truncation error converges more slowly 
than the discretisation error. For example, it may be the case that G 0(h 2 ) while G 0(h) or even 
Th G 0(1) meaning that the truncation error does not reduce with grid refinement. This has been observed 
by several researchers, e.g. [23-25, 11, 26-28]. Unfortunately, such large truncation errors appear not 
only on irregular grids but also on smooth grids under the usual discretisation schemes near the domain 
boundaries [25] and near the interfaces between different levels on locally refined grids [11]. In all these 
cases, when second order accurate discretisation schemes are used, numerical experiments show that the 
high truncation errors of order 0(1) do not cause degradation of the convergence rate of the discretisation 
error which remains 0(h 2 ). An interesting possible explanation for this based on arguments from potential 
theory is sketched in [29]. More rigorous mathematical arguments can be found for the domain boundaries 
e.g. in [30] and [31]. This phenomenon requires some investigation in relation to the use of the truncation 
error as a local refinement criterion, because in the regions of high truncation error it seems likely that 
there will result perpetual refinement without actual benefit in terms of reduction of the discretisation 
error. 

3. Truncation error estimate 

A popular method for estimating the truncation error works by transferring (restricting) the solution 
obtained on a given grid h to a coarser but geometrically similar grid i7, and applying the discrete operator 
of the coarse grid (constructed using the exact same discretisation schemes as the operator of the fine 
grid h) to this restricted solution [4, 32, 8, 33, 10, 25]. This concept arises naturally in FAS multigrid 
methods (Full Approximation Storage schemes - see any of [4, 34, 8]), from which the estimate originates. 
The method is valid under the assumption that the truncation and discretisation errors vary smoothly 
within the domain, and there exist continuous functions r and e such that: 


r h = i%[hp-r + om\ ( 3 ) 

e h = 7* [h? ■ e + 0(h q )\ (4) 

where p is the order of the finite volume method and q > p, with both p and q being integers. Iq is an 
operator that discretises a continuous function by keeping only the values of the function at the centres 
of the CVs of grid h. 

Suppose as discussed in Section 2 that to solve the differential equation N(u) = b the problem is 
discretised on a grid h using the finite volume method to obtain the algebraic system Nh(uh ) = bh by 
neglecting the truncation error, which is solved to obtain an approximate solution u^. Suppose also that 
there is available a grid 2 h which is geometrically similar to grid h but which has everywhere twice the 
grid spacing. Consider the following systems on grid 2/i, whose solutions are increasingly more accurate: 


N2h(u2h) = hh (5) 

N 2h{u2h) = hh ~ A = N 2 h {u2h) ~ A ( 6 ) 

N 2 h(u 2 h) = hh - T 2 h = N 2 h(u 2 h) - T 2 h (7) 

where 

A = hh - N 2h ( I 2 h h u h ) (8) 

is called the relative truncation error on grid 2 h with respect to grid h. In (8) uk is the approximate 
solution obtained on grid h by the finite volume method, and iff 1 is a transfer (restriction) operator which 
transfers the function u h to grid 2 h. For the time being it will be assumed that this transfer does not 
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introduce any errors, in the sense that, for example, iffffu^ = U 2 h (uh = IqU and U 2 h — iff^u being the 
values of the exact solution u of the differential equation at the centres of the CVs of grids h and 2 h 
respectively). This is a natural assumption if the CV centres of grid 2 h are a subset of the CV centres 
of grid /i, as in node-centred grids, and direct injection is used (i.e. the values of U 2 h are simply copies 
of the corresponding values of uh at the coincident points). But in cell-centred grids like those used in 
the present study this assumption is not valid, and the implications will be discussed later on. The finite 
volume operator N2h is constructed on grid 2 h using the exact same discretisation schemes as on grid 
h. 

The following can be noted about the solutions of systems (5) - (7): U2h is the approximate solution 
on grid 2 h obtained by solving (5). By substituting (8) into (6) one can see that the solution of (6) is 
u 2h — iffffuh, i.e. the solution of grid h restricted to grid 2 h. So, by subtracting r 2h from the right hand 
side of (5) one gains the accuracy of grid h on grid 2 h. Finally, by subtracting the actual truncation error 
T 2 h from the right hand side of (5) one obtains the exact solution of the PDE ( U 2 h ) by solving (7). 

Let N 2h be the Jacobian matrix of evaluated at U2h (N^h = N 2 h if ^2h is linear), and assume 
that the functions U2h, and u 2h are close enough such that the following hold: 

N 2h (if %) = N 2h (u 2h ) + N' 2h ■ (if % - u 2h ) + O (5?) (9) 

^2h(u2h) = N 2 h{u2h) + A^2 h * ( u 2h ^2h) + O (8ff) (10) 

where Si — iffffun — U2h and S 2 — U 2 h ~ &2h- Now, Si and S 2 are of the order of the discretisation error 

(see also equation (13) below) and so 0(S\),0(S 2 ) G 0(h 2p ) due to (4). By rearranging (9) - (10) and 

comparing to (6) - (7) respectively one obtains that: 


N' 2 h -(l 2 h h u h -u 2h ) = —r 2h + 0(h 2p ) 

N 2 h ■ ( u 2 h ~ U 2h ) = —T 2h + 0 (h 2p ) 


Now, u 2 h - u 2h = e 2 h, and also: 


( 11 ) 

( 12 ) 


if % ~ U 2 h = ( U 2h - U 2h ) - (u 2h - Ih k Uhj = t2h ~ Ih he h 
= if [(2 h) p e + 0(h q )\ - iff [h p e + 0{h q )] 

= I$ h [2 p h p e-h p e + 0{h q )} 

= (l-L^)/f [2 p h p e + 0{h q )] 

= (^r) ^ + om 


(13) 


where (4) has been used for both and 62 h, an d also it has been assumed that iff 1 has not introduced 
any errors (iff 1 1 q = iff 1 ), as stated above. Therefore, (11) and (12) become: 


+ N ^-°( hq ) = -4 + 0{h 2p ) (14) 

N' 2h ■ e 2h = —T 2h + 0(h 2p ) (15) 

Usually it will be the case that 2 p > q. Also, it may be assumed that N 2h • 0(h q ) G 0(h q ) because neither 
N 2h nor the vector 0(h q ) change arbitrarily as the grid spacing, /i, decreases, but they approximate 
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specific continuous counterparts. The eigenvalues of N 2h are therefore expected to not change much with 
grid refinement. So, substituting (15) into (14) and rearranging one gets: 

9 p 

T 2h = j ~ T 2h + °( hq ) ( 16 ) 

Finally, (3) implies that: 

T h = F I$ h T 2h + 0(h q ) (17) 

where I 2h is a transfer (prolongation) operator which transfers a function from grid 2 h to grid h. This, 
combined with (16) and the definition of r 2h (8) gives the final result: 


Th 


1 jh 

2p — 1 2h 


b 2h ~ N 2h ( P h h u h 


.)] + 0(h q ) 


_l n(M\ 


(18) 


The above estimate differs from the actual truncation error by 0(h q ). Since the truncation error 
itself has a magnitude of 0(/i p ), the relative error in the approximation is (r^ — f^)/r^ G 0(h q ~ p ). In 
most cases it will be the case that q — p — 1, but there are some central difference approximations where 
due to cancellation of every second term in the corresponding Taylor series it happens that q — p = 2. So, 
it may be stated that the approximation fh (18) is either first- or second-order accurate. 

The analysis so far has assumed that the transfer operators iff 1 and I 2h are “perfect” in the sense that 
they introduce no errors. In practice this is not possible on cell-centred grids, and so an attempt will now 
be made to examine the effect of the order of accuracy of these operators on the order of accuracy of the 
estimate (18). Starting with the prolongation operator, suppose an r-th order accurate operator I 2h such 
that = I h 2h u 2h + 0(h r ), / 2 \ being again the “perfect” operator which introduces no errors. Then 

I% h T2h = l2h T * h + 0(h p+r ) because r 2h G 0(h p ). So (17) becomes r h = ^ h r 2h + 0{h p+r ) + 0(h q ). This 
means that the order of accuracy of the estimate (18) will degrade only if p + r < q, or r < q — p. If 
q — p — 1 then using a first-order accurate operator (r = 1) suffices. In the present study a first-order 
accurate prolongation operator is used for this purpose which sets the value at each CV of grid h equal 
to the value of its parent CV of grid 2 h. This causes the estimate (18) to predict the same truncation 
error for all sibling CVs, i.e. CVs which share the same parent. In general therefore all siblings are 
refined together. But this is already a requirement of the refinement procedure in order to ensure that 
the underlying grid 2 h always exists (see Section 4). 

Now, suppose an r-th order accurate restriction operator I^ h such that I^ h Uh = Iy^ u h + 0(h r ). Then 
(18) can be written as: 


Th = 


1 jh 

2p — 1 2h 

1 rh 

2p — 1 2h 

1 jh 

2p — 1 2h 


P h h u h + 0(h r ))] + 0(h q ) 


= n + 0(h min 


b 2 h ~ N 2 h (. 

b 2h - N 2h ( I 2 h h u h ) - iV' fe • 0(h r ) + 0{h 2r )J + 0(h q ) 
b 2h - N 2h (j P h h u h )] + 0(h r ) + 0{h q ) 


(19) 


So now ( Th — Vi)/rh G 0{h mn( T^~ p ) so that in order for the estimate fh to converge to the actual it 
is necessary that r > p. The maximum accuracy is obtained for r > q. These theoretical findings are 
demonstrated in numerical experiments in [25]. For the present study, the better of the two restriction 
operators tested in [25] is used, which is third-order accurate on smooth grids and second order accurate on 
irregular grids (including near domain boundaries and near interfaces between different local refinement 
levels). 
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In the case of non-smooth grids, the accuracy of the restriction operator is not the only issue of 
concern. Indeed, in this case (3) does not hold (although (4) does, in general, hold) so (18) does not 
hold because p is wrong or varies across the domain. Yet even in this case (18) should be able to roughly 
predict the order of magnitude of the truncation error, due to the validity of (4), which could be sufficient 
for local refinement purposes. In this case one can regard (18) as calculating the truncation error on grid 
2 h by applying the coarse grid operator N^h to the fine-grid solution instead of the exact solution which 
would give the exact truncation error. The fine grid solution may be regarded as a good approximation 
to the exact solution on grid 2 h due to (4). 

Finally we note that the estimate (18) can also be used with any other coarse-to-fme grid ratio, 
H = p • /i, by substituting p instead of 2 in the denominator. 


4. Equations and discretisation 


For simplicity, the 2-dimensional steady-state incompressible Navier - Stokes equations with constant 
density and viscosity are solved in this work. Integrated over a CV P, the equations solved are: 
x-momentum: 



y-momentum: 



continuity: 



( 20 ) 

( 21 ) 

( 22 ) 


where Sp is the boundary surface of CV P, dS is an infinitesimal element of this area pointing out of the 
CV, V = (u,v) is the fluid velocity, p, p, and p are the density, viscosity and pressure respectively, and i, 
j are the unit vectors in the x and y directions respectively. The second terms on the right-hand sides of 
the momentum equations are zero for constant density and viscosity, yet they have been retained because 
their discrete counterparts were present in the Fortran code that was used for this study (in their discrete 
form these terms are not zero but they add to the truncation error, tending to zero as the grid is refined; 
cf. end of Sec. 5.3). 

The in-house code used for the present study can handle grids composed of quadrilateral CVs. A 
structured or block-structured grid is used as a starting point, and local refinement levels are added 
on top of that by the adaptive procedure. Refinement of a CV takes place by splitting it along the line 
segments which join its centre with the centres of its faces, into four child CVs. The original CV, called the 
parent CV, is retained in the data structure, and is used by a multigrid algorithm to accelerate algebraic 
convergence. Each CV belongs to a certain grid level : these levels are numbered, and the children of a 
CV of level l belong to level l + 1. Each level consists of global and local CVs: Global are those which 
do not have children, and local are those which do have (see Figure 1). The actual grid on which the 
problem is discretised and solved at any one time (composite grid) is the set of all global CVs of all levels. 
The CVs of the composite grid, although geometrically quadrilateral, are treated as logically polygonal 
because one or more sides of a CV may actually consist of two faces, if the CV borders a finer level on 
that side. A face is defined as the piece of the boundary of a CV which separates it from another single 
CV. For example, the left side of the CV whose centre is N in Figure 2 consists of two faces, one of which 
is face / separating it from the CV whose centre is P, which belongs to a finer level. Neighbouring CVs 
are not allowed to be more than one level apart. For more information on the grid structure see [11]. 

The discretisation schemes are those used in [25] and are similar to those described in [22]. They are 
central difference schemes with correction terms for grid non-orthogonality, skewness etc. Suppose a face 
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Grid level 1 level 2 level 3 level 4 

Figure 1: A composite grid and its decomposition into levels. The local CVs of each level are marked in grey. 



/ of area Sf and centre c separating two CVs whose centres have position vectors P and N as shown in 
Figure 2. Suppose also c f is the point on line PN which is closest to c. Also let points P', N' be such 
that the line segment P'N' is perpendicular to /, its length is equal to the length of PN, and its centre 
is at c. The convective and viscous terms of (20) are discretised thus: 


where: 


1 puV • dS ~ F t c 

(23) 

Js f 


J^nVu-dS K nSf 

(24) 

Uh,c = (Uh)d + (VftUfcV • (c - d) 

(25) 

u h,P' = u h,P + (Vh u h)p ■ (■ P' ~ P) 

(26) 


In the above equations, Uh is the grid function of the values of u at the CV centres of grid h and Uh,p is 
the value of that function at the centre of CV P. The overline denotes linear interpolation at point c' from 
the values at points P and N. The discrete gradient operator is based on a weighted least squares 
method that is second-order accurate on smooth grids and first-order accurate on irregular grids (such as 
the one shown in Figure 2), described in [13, 35]. The mass flux F^j through face / is calculated using 
a variant of momentum interpolation to avoid pressure oscillations which may otherwise be exhibited by 
the solution when both velocity and pressure are stored at CV centres [22]: 


PS} 


f 


Fhj = pVh, c ■ Sf + ( Ph,p ~ Ph,N ) - (VhPh) p±n ■ (P - N) 


(27) 


Here the overline denotes interpolation of the gradient at the middle of the line segment PN. The 
original concept of momentum interpolation was proposed in [36], but it has the shortcoming that Af is 
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Figure 3: A boundary CV with its boundary face b. 


calculated in the context of the SIMPLE solution method [37]. To separate completely the discretisation 
scheme from the algebraic solution method, a modified momentum interpolation was proposed in [25], 
where Af is calculated directly from flow variables and grid geometry: 


A 


pSf\Vh,c m nf\ + pSd\Vh, c -ROT(ftf)\ + 2/i 


~ S 1, S A 

U, S f \ 


(28) 


Here tif is the unit vector normal to face /, ROT is a function which rotates a vector by 90 degrees, and 
Sd = (N — P) • ftf. The pressure term of (27) does not allow the development of pressure oscillations in 
the solution. It is an artificial term totally belonging to the truncation error, and, as discussed in [25], its 
contribution to the truncation error is 0(/i 4 ) on smooth grids and 0(/i 3 ) on irregular grids. Therefore it 
does not affect the overall accuracy of the method. 

The convection approximations (23) of all faces of a CV contribute 0(h ) to the truncation error 
on non-smooth grids, and 0{h 2 ) on smooth grids because in the latter case parts of the contributions 
cancel out between opposite faces of the CV, as shown e.g. in [35]. Also shown in [35] is that the 
diffusion approximations (24) contribute 0(1) to the truncation error on non-smooth grids and 0(/i 2 ) on 
smooth grids. Nevertheless, as discussed in Section 2 , the discretisation error remains 0(h 2 ) under these 
discretisation schemes, even when the truncation error is 0{h ) or 0 ( 1 ). 

The discretisation of the pressure terms of (20) - ( 21 ) is similar to that of the convection terms. As 
for the second viscous terms (second terms on the right hand sides of ( 20 ) and ( 21 )), they are discretised 
using the midpoint rule approximation, using the value of the gradients at point d as obtained by linear 
interpolation. 

It is appropriate also to discuss briefly the implementation of the wall boundary condition, which is 
the only boundary condition used in the test case of Section 6 . The original version of the CFD code used 
for the present study followed the popular approach described in [ 22 ] where instead of just implementing a 
Dirichlet boundary condition for the momentum equations, the continuity equation is also used to further 
refine the boundary condition. For example, Figure 3 shows a boundary CV with a boundary face / 
oriented along the x direction. If the wall is rigid and moves with a constant tangential velocity then 
du/dx = 0 implies that also dv/dy = 0. Therefore the normal viscous stress is zero at a wall, and so only 
a tangential viscous force is implemented on the boundary, discretised as f Sb —fi^dS ~ • 

But if the tangential velocity of the wall is not constant (imagine the wall as being elastic, stretching 
at some parts and contracting at others) then du/dx 7 ^ 0 implies that also dv/dy 7 ^ 0 and there is also 
a normal viscous stress. To allow for this the CFD code was modified to include also a normal viscous 
force: f Sb —fi^dS ~ ~where vb — 0 on the wall (the normal component of velocity is zero, 
but its normal derivative is not, in general, zero). The 2 nd viscous terms on the right hand sides of (20) 
- ( 21 ) are also discretised by taking the gradients ( X7u)b and (Vv)b as equal to the respective gradients 
at the centre of CV P. Preliminary numerical experiments on the classic lid-driven cavity problem with 
rigid walls, for which both approaches are valid, showed very little difference in the results (the method 
of [22] produced very slightly better results, as evaluated against the benchmark results of [38]). 

The discretisation of (20) - ( 22 ) for each CV results in a system of algebraic equations whose unknowns 
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are the velocity and pressure at each CV centre. This system is solved using the popular SIMPLE 
algorithm [37] in a multigrid context which is such that it allows for the existence of local refinement 
levels that do not extend throughout the domain. This multigrid method was proposed in [11]. SIMPLE 
is a fixed-point iteration type algorithm, where during each outer iteration a linear system is solved for 
each of the velocity components, followed by a linear system for the pressure. The velocity component 
linear systems are obtained by linearising the corresponding momentum equations using values from the 
previous SIMPLE iteration in the nonlinear terms, and in terms containing variables other than that 
particular velocity component. The velocity field that results from the solution of these linear systems 
does not satisfy the continuity equation in general, and so a Poisson-type linear system is constructed 
and solved to calculate corrections to the pressure field that would force the velocity field to be more 
mass-conserving. This completes one outer SIMPLE iteration, and the whole procedure is repeated in 
the next outer iteration. This brief description of the SIMPLE algorithm will be useful to the unfamiliar 
reader as mention of it is made a few times in this paper, although the methodology and results are in 
no way connected to a particular algebraic solver. 

The SIMPLE-multigrid procedure solves the algebraic system up to a residual. To examine the effect 
of this on the accuracy of the solution, consider again the abstract problem N(u) = b. As already noted, 
the exact solution uh of the discrete system which approximates the continuous problem satisfies equation 
(2), repeated here: 


Nh{uh ) = N h (u h ) + r h (29) 

where uk is the exact solution of the continuous problem. Now, the iterative procedure does not solve 
for Uh exactly, but suppose that after the k -th iteration the algebraic residual is and the approximate 
solution is u\: 


N h (u k h ) = N h (u h ) + r k h = N h (u h ) + (r h + r k ) (30) 

where (29) has been used to obtain the second equality. Comparing with (29), this result shows that 
Ufr corresponds to a “truncation error” + r^, just like uh corresponds to the truncation error r^. To 
obtain the maximum accuracy, iterations must be carried on until the residual has dropped below the 
truncation error r^. This is important also in order to accurately estimate the truncation error using 
(IS) 1 . 

5. Local refinement 

5.1 . General procedure 

In the present work local refinement takes place as follows: After solving the equations on a given 
grid (which may be composite), the truncation error is calculated at every CV using equation (18). This 
truncation error estimate is used to calculate another quantity which will be used to select the CVs to be 
refined (such quantities will be described shortly). Then all CVs are sorted in descending order according 
to the absolute value of this quantity. A fraction R of the CVs, those with the largest value of this 
quantity, are marked for refinement, i.e. the first R • K CVs of this descending list are marked, where K 
is the total number of CVs in the grid. An efficient algorithm should be chosen to perform the sorting 
because an inefficient algorithm could require a significant amount of time if the grid is large in terms of 
number of CVs. This has been observed in the present work when the selection sort algorithm was first 
used due to its simplicity. The situation improved greatly when the quicksort algorithm was used instead, 


1 A detailed discussion of the effect of the iteration error on the truncation error estimate is given in [39], where it is shown 
that it is possible, using specialised restriction operators, to obtain a good estimate of the truncation error even before the 
algebraic residuals have dropped to a low level. The results of [39] were not known to the authors when the present study 
was conducted. 
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Figure 4: Local Refinement, (a): Suppose a single CV (•) is marked for refinement according to the criterion, (b): its 
neighbours (o) are also marked, as a safety margin, (c): The siblings (■) of any marked CVs are also marked, as required 
for the existence of the underlying grid, (d): After all CVs have been marked, refinement takes place to produce the new 
grid h, whose underlying grid 2 h is (e). 


with the cost of sorting becoming insignificant compared to the overall computational cost (see e.g. [40] 
for a description of the algorithms). 

Marking of the CVs takes place as shown in Figure 4: When a CV is marked, its neighbours are also 
marked as a safety margin (Figure 4 (b)). Also, as shown in Figure 4 (c), for each CV marked its siblings 
are also marked (i.e. all CVs which share the same parent CV), which is necessary so that the new grid h 
which will result when refinement takes place will have an underlying grid 2 h on which the new truncation 
error will be calculated. This marking procedure results in more than R • K CVs being marked overall. 
Furthermore, if a system of PDEs is solved, like in the present case where there are 3 equations for each 
CV (2 momentum equations (20) - (21) and the continuity equation (22)), the marking procedure can be 
repeated for as many of these equations as desired, and the results are logically OR-ed to establish the 
CVs that are to be refined. 

After all relevant CVs have been marked, refinement takes place by splitting each marked CV into 
four children, at the line segments joining their centres with the centres of their four sides, as shown in 
Figure 4 (d). By marking all siblings together, this procedures ensures the existence of the coarse grid 2 h 
(Figure 4 (e)) which consists of the parents of all the global CVs, or equivalently of all the local CVs whose 
children are global. The equations are then discretised and solved on the new grid /i, the truncation error 
is calculated using the coarse grid 2/i, and the refinement procedure is repeated to produce yet another 
grid. This cycle is repeated a fixed number of times. 

5.2. Quantities used as refinement criteria 

In the present work, the quantities that are used to assess whether a CV P should be refined or not 
are the absolute values of the following: 

1. The truncation error itself (t/^p). 

2. The truncation error times the volume of the CV (t/^p • 5ftp). 

3. The truncation error divided by the corresponding main diagonal coefficient of a matrix A^ which 
comes from linearising the discrete operator N^ (r^p/A/p^p). 

The quantity 1 is used for example in [7, 9, 10]. The quantity 2 is used e.g. in [8]. And quantities 
similar to 3 are used e.g. in [13, 14] (except that there the residual is used instead of the truncation error). 

The reasoning behind quantity 3 can be sketched as follows: If Nk is linear, then the relationship 
(2) between the truncation and discretisation errors simplifies to N^e^ = — Tp. Then, dividing r h by the 
main diagonal coefficients of N\ is equivalent to performing a single Jacobi iteration on that system with 
zero initial guess. This produces a very rough estimate of the discretisation error e^ produced by the 
truncation error in each CV that takes into account the physical processes (e.g. convection, diffusion) 
inherent in N^. This rough estimate of e^ is actually used as the refinement criterion. 
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The main diagonal coefficients used in quantity 3 are those of the velocity linear systems in the SIMPLE 
solution method, as long as the systems are expressed per unit volume and the matrices of coefficients 
are constructed using the UDS (first order upwind) scheme for convection with any higher-order schemes 
being implemented through deferred correction [41, 22]. The UDS ensures that only the flow out of the 
CV is taken into account, and not the flow in that would lead to cancellations and possibly near zero or 
negative contribution of convection to the coefficient [13, 14]. 

If the systems are expressed per unit volume then quantity 3 for a CV P becomes |r^p/App| and 
\T][p/A v pp\ for the x- and y- momentum equations respectively, A u , A v being the matrices of coefficients 
of the linear systems for the velocity components. If the systems are not expressed per unit volume then 
quantity 3 is calculated as \(r% p 5Qp)/Ap p \ and \ (r% p 5Qp)/A pp \. The coefficients are readily available if 
the SIMPLE solution method is used, otherwise they have to be calculated. When not expressed per unit 
volume their form is similar to (28), only that (28) is staggered, i.e. centred around a face centre, while 
4 p ; p would be centred around the centre of CV P. An examination of (28) shows that grid refinement 
causes Ap^p to tend to a constant; its convection part tends to zero because it is proportional to the 
CV dimensions Sf and S^, while the viscous part is constant because the refinement procedure adopted 
here does not alter the aspect ratio Sf/Sd of the CVs. If all the CVs of the grid have the same aspect 
ratio then as the grid becomes finer Ap^p tends to the same value for all CVs (for constant y) and the 
quantities 2 (r/^pADp) and 3 (rh^p-SQp/Ap^p) become equivalent. If the CV aspect ratio varies across 
the grid then criterion 3 has a preference for refining CVs with an aspect ratio close to 1. 

In the present work quantity 3 is not calculated for the continuity equation because a corresponding 
coefficient is not readily available, although it should be possible to derive an approximate equation 
relating the truncation error of the continuity equation to the discretisation error of pressure in the same 
way that the “pressure correction equation” is derived in the SIMPLE method. 

5.3. Internal boundary treatment 

There is another important issue which must be addressed. The derivation of the truncation error 
estimate in Section 3 has been based on a number of idealised assumptions, not all of which are valid 
in a realistic case. Specifically, as discussed in Section 2, the truncation error becomes 0(1) at domain 
boundaries and at the interior boundaries between different grid levels. This is because although the 
grid is generally Cartesian, it exhibits all kinds of distortions directly on the level interface as shown in 
Figure 2. Therefore the truncation error distribution is discontinuous at these regions, and the index p 
of equation (3) is zero there, while at the smooth interior regions of the grid it equals 2. Yet, as also 
discussed in Section 2, equation (4) remains valid even on non-smooth grids. Even though (3) is not 
completely valid on composite grids, equation (18) should provide a rough estimate of the truncation 
error, due to the validity of (4). That is, since the truncation error estimate works by assuming that the 
solution of the fine grid h is roughly equivalent to the exact solution of the PDE when considered with 
respect to grid 2/i, and since this assumption is roughly true even on distorted grids due to (4), should 
provide a rough estimate even at distorted regions of the grid, albeit not as accurate as at smooth regions 
where all the assumptions hold. 

Figure 5 shows an example of the high truncation errors that occur near grid level boundaries. The 
x-momentum equation truncation error estimate is displayed along with the underlying grid 2 h for one of 
the test cases which will be presented in Section 7. It is a good idea to display the coarse grid 2 h rather 
than grid h because most of the truncation error calculation occurs on grid 2/i, and since an injection 
prolongation operator is used all children of a CV of grid 2 h inherit the same value of truncation error 
calculated at their parent. It may be seen that the truncation error is several orders of magnitude higher 
near the level interfaces than in the interior of the levels. This is expected since in the former case 
Th G 0(1) while in the latter case G 0(h 2 ). It may also be noticed that the high truncation error 
regions appear to extend up to two CVs of grid 2 h away from the level interface, i.e. up to 4 CVs of grid 
h. This is an artificial result of the truncation error estimate (18) which calculates the truncation error 
mostly on grid 2/i, caused by the discontinuity of the truncation error distribution. In fact the actual high 
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Figure 5: Truncation error (estimate) of the ^-momentum equation at a region with grid level interfaces. The underlying 
grid 2 h is shown. 


truncation error occurs only up to two CVs of grid h away from the interface. A similar complication 
also seen in the Figure is that in the high truncation error regions the truncation error often appears to 
alternate in sign between neighbouring CVs of grid 2 h. This should be expected to occur in fact between 
neighbouring CVs of grid /i, but a usual prolongation operator cannot produce this. 

Therefore, to transfer correctly the truncation error estimate from grid 2 h to grid h a specialised 
prolongation operator should be used that takes into account the similarity between CVs of grid 2 h and 
CVs of grid h. For example, the truncation error at CVs of grid h that are in the second row away from 
the level interface should be based not on the calculations performed at their parents (which are in direct 
contact with the level interface) but on their parents’ neighbours that also lie on the second row (of grid 
2 h) of CVs away from the interface. Anyway, in the present work such a specialised prolongation operator 
is not used, because the issue is not how to transfer accurately the 0(1) truncation error from grid 2 h 
to grid /i, but how to deal with it in the context of local refinement, which assumes that refinement 
of CVs causes a reduction in truncation error. In this case, of course, such a reduction does not occur 
because E 0(1), so refinement appears to bring no benefit, while on the other hand it increases the 
computational cost by increasing the number of CVs. Three different strategies to deal with this issue 
are tested in the present work: 

1 . Take no particular action. CVs near grid level interfaces may be refined just as all other CVs of the 
grid. 

2 . Do not allow refinement of CVs that are near level interfaces. This includes CVs that are up to 4 
positions away from the interface, as discussed above. 

3. Do not allow refinement of CVs that are on the fine side of the interface, but allow refinement of 
CVs that are on the coarse side of the interface. This does not allow the creation of an even finer 
level, but rather the extension of the existing fine level towards the coarse side. 

These strategies will be tested in Section 7. Similar problems occur at the computational domain 
boundaries where the boundary conditions are applied, but in the present work no special restrictions on 
CV refinement have been applied there. 

To check whether a CV is within a 4-CV width zone on either side of the level interface, it is checked 
whether its parent or any of the parent’s neigbours of the same level as the parent are in contact with the 
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level interface (the 4-CV width zone appears as a 2-CV width zone on the underlying grid). The parent 
is checked first; since it belongs to the underlying grid it will be in touch with the interface if any of its 
neighbours of the same level do not belong to the underlying grid, i.e. if they either do not have children, 
or they have grandchildren. If all the parent’s neighbours belong to the underlying grid then the parent 
is not on the interface, but the check must be repeated for its neighbours as well. Alternatively, one could 
traverse the linked list of faces which lie on the interface (they separate a fine from a coarser CV). The 
data structures are implemented using pointers in Fortran 95. 

Closing the present section, we would like to make a couple of comments. The first concerns the reason 
why the high truncation error extends also to the second CV away from the interface. The reason is that, 
according to formulae (23) - (26), the calculation of the fluxes through a face that separates such a CV 
from another CV that touches directly the level interface involves also the gradients at this neighbour 
CV, which are only first order accurate, as mentioned in Section 4 (because the grid is distorted directly 
on the interface, as shown in Figure 2 ). On the other hand, one may notice that the gradients in the 
terms (23) - (26) are used only if the grid exhibits skewness, non-orthogonality or unequal spacing of the 
CV centres from that face. However, even in the absence of these geometric distortions the gradient at 
the CV which touches the interface is still used to calculate the second viscosity terms of (20) - ( 21 ), i.e. 
the second terms on the right-hand sides, for the face that separates the two CVs. So even if the grid is 
smooth, a CV that has a neighbour in direct contact with the level interface will exhibit large truncation 
errors. 

The second comment concerns the alternating sign of the truncation error at neighbouring CVs near 
the level interfaces. This is due to the fact that when a conservative discretisation scheme is used, like in 
the present work, each face contributes equal and opposite contributions to the equations of the CVs on 
either side of the face. This holds also for the truncation error contributions: the discretisation of the flux 
through a face contributes equal and opposite contributions to the truncation errors of the CVs on either 
side. Therefore the total truncation error in the domain is determined only by boundary faces and source 
terms. Also, if a pair of neighbouring CVs has 0(1) truncation errors of similar magnitude but opposite 
sign, then the effect of these truncation errors (in terms of the discretisation errors they produce) at distant 
locations will tend to cancel out. Local refinement of these CVs, although it will not reduce the truncation 
error since it is of 0(1), will produce smaller CVs of opposite truncation error that are closer together. 
Bringing the two sources of opposite strength closer together would cause the overall discretisation error 
that they produce at any given location to decrease. So local refinement could prove beneficial even in 
cases when E 0(1), which is in line with the observation mentioned in Section 2 concerning the validity 
of equation (4) even in cases when (3) does not hold. For example, in [11] a numerical experiment is 
described where the lid-driven cavity problem is solved on a series of similar composite grids of increasing 
fineness (each successive grid is produced by refining every CV of the previous grid) and although the 
truncation error does not decrease near level interfaces, the discretisation error does decrease as 0{h 2 ). 

6. Description of the test case and of the evaluation methodology 

6.1. The regularised lid-driven cavity case 

The refinement strategies previously presented were tested on a lid-driven cavity problem. The classic 
lid-driven cavity problem is the most popular problem for testing and validating new techniques for 
incompressible flows (see e.g. [42, 38, 43, 11] and references therein). It concerns the flow in a square 
cavity with three stationary sides and one side (lid) - usually taken as the top side as shown in Figure 
6 - moving at a constant tangential velocity which drives the flow. Unfortunately the resulting flow 
field is discontinuous because the velocity changes abruptly from a non-zero constant value (top lid) to 
zero (side walls) at the top corners of the cavity. This causes very large truncation errors at these two 
corners, which increase with refinement (see e.g. [11]). This would pose an additional complexity for the 
present investigation of the efficiency of truncation error-based local refinement, making it more difficult 
to interpret the results. So, to keep things simple, the so-called regularized lid-driven cavity problem 
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U = - 16x 2 (l-x) 2 
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Figure 6: The 32 x 32 grid and the location of the probe points for the lid driven cavity problem. 

[44, 45] was selected instead, where the tangential velocity of the lid is not constant (the lid may be 
viewed as made of an elastic material that can stretch and contract), but is given by a polynomial that 
is such that the velocity and its first derivative are zero at the edges of the lid: 

U{x) = — 16x 2 (l — x) 2 (31) 

It is necessary that dU/dx be zero at the lid edges so that the continuity equation is satisfied: At the 
lid corners it holds that dv/dy = 0 because v = 0 at the side walls, so it must also hold that du/dx = 
dU/dx — 0. Equation (31) sets the maximum velocity at the centre of the lid with \U max \ = 1 m/s. 
There are no Dirichlet boundary conditions for pressure but instead pressure is linearly extrapolated to 
the boundaries from the interior [22]. After solution, the whole pressure field is adjusted by a constant 
so that p = 0 at the centre of the cavity. 

First, the problem was solved on a series of structured grids with uniform grid line spacing: 32 x 32 
CVs (shown in Figure 6), 64 x 64 CVs, etc. up to 2048 x 2048 CVs. These will be referred to as structured 
grids in the following, while the locally refined grids will also be referred to as composite grids. Then 
using Richardson extrapolation (see e.g. [22]) a very accurate solution was obtained which was treated as 
the exact solution. This was repeated for a range of Reynolds numbers: Re — 100, 1000 and 10000 where 
Re — pUrn ^ xL , The values of the constants are p — 1 kg/m 3 (density), Umax = 1 m/s (max. lid velocity) 
and L = 1 m (length of cavity side). The viscosity p was varied to obtain the desired Reynolds number. 

For the classic lid-driven cavity problem, published results (e.g. [43]) indicate that the flow is steady 
for Re = 100 and 1000, but not for 10000 where it becomes periodic in time. For the regularised lid- 
driven cavity problem it is suggested in [44, 45] that the flow is steady at Re — 10000, although it becomes 
unsteady at slightly higher Reynolds numbers. The present work deals only with steady-state flows, so the 
flow at Re = 10000 is considered to be steady (which is most probably the case as suggested in [44, 45]). 

Multigrid cycles were carried out on each grid until the residuals of all equations in every CV had 
dropped below 10 -8 to ensure that they were well below the truncation error (see eqn. (30)) 2 . For 
Re— 100, iterations converged on each grid (including composite grids) in about 13 V(2,2) cycles. For 
Re — 1000, converge was obtained with around 12 W(2,2) cycles (more on coarse grids). And for Re — 
10000 using again W(2,2) cycles, convergence was obtained with 12 and 10 cycles for the 1024x1024 and 


2 Actually it suffices to ensure that the residuals have dropped below the truncation error by, say, an order of magnitude, 
according to equation (30) . Any further reduction would be a waste of computational effort because the accuracy would not 
increase. In this respect the 10 -8 threshold that we used throughout the domain is supererogatory. 
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2048 x2048 grids, respectively, while for coarser grids more cycles were necessary. Finally we note that 
for i?e=1000 and i?e=10000 within each cycle the corrections were smoothed before being prolonged to 
the finer grid according to the procedure described in [11]. 

Tables 1 and 2 present the “exact” solution as obtained using Richardson extrapolation at a number of 
points along the vertical and horizontal centrelines respectively. The points are marked on Figure 6, and 
they are the same points used in numerous studies on the classic lid-driven cavity problem; they were first 
used in [42], and also later by other authors such as [38, 43] who compared their results against those of 
[42]. So we present here the respective results for the regularised cavity case, to be available as benchmark 
results for future studies. The present results were obtained by first linearly interpolating the solution 
of the 1024x1024 and 2048x2048 grids at these points and then performing Richardson extrapolation, 
assuming 2 nd order convergence. Repeating the procedure using grids 512x512 and 1024x1024 instead 
showed no difference up to the 7 digits shown in Tables 1 and 2 for Re=100 and Re=1000 except for a 
couple of points near the lid. For Re=10000 the last digit (or the last couple of digits near the lid) may 
not be accurate. 


Table 1: Numerical solutions obtained along the vertical centreline. 


y 

Re = 

u 

100 

p 

Re = 

u 

1000 

V 

Re = 

u 

10000 

V 

0.0000 

0.0000000 

0.0250098 

0.0000000 

0.0607925 

0.0000000 

0.0605268 

0.0547 

0.0306724 

0.0250583 

0.1054923 

0.0605163 

0.3178112 

0.0557021 

0.0625 

0.0345378 

0.0250546 

0.1184943 

0.0603569 

0.3252990 

0.0539831 

0.0703 

0.0383053 

0.0250475 

0.1312872 

0.0601478 

0.3243764 

0.0521744 

0.1016 

0.0525963 

0.0249789 

0.1813140 

0.0586509 

0.2920120 

0.0451693 

0.1719 

0.0815204 

0.0244648 

0.2668322 

0.0496716 

0.2427313 

0.0316099 

0.2813 

0.1222629 

0.0217104 

0.2303403 

0.0257505 

0.1690492 

0.0151767 

0.4531 

0.1630748 

0.0070318 

0.0873616 

0.0028316 

0.0529534 

0.0012714 

0.5000 

0.1612522 

0.0000000 

0.0519093 

0.0000000 

0.0209067 

0.0000000 

0.6172 

0.1165433 

- 0.0211446 

- 0.0397520 

- 0.0010470 

- 0.0612785 

0.0014953 

0.7344 

0.0140787 

- 0.0389617 

- 0.1399778 

0.0063314 

- 0.1496488 

0.0092355 

0.8516 

- 0.1591120 

- 0.0426021 

- 0.2387008 

0.0201424 

- 0.2503312 

0.0216989 

0.9531 

- 0.5752639 

- 0.0250529 

- 0.3131419 

0.0304935 

- 0.3470263 

0.0327463 

0.9609 

- 0.6336426 

- 0.0223860 

- 0.3464474 

0.0310851 

- 0.3473972 

0.0333507 

0.9688 

- 0.6980380 

- 0.0194587 

- 0.4019992 

0.0317165 

- 0.3435578 

0.0338717 

0.9766 

- 0.7668059 

- 0.0163600 

- 0.4880689 

0.0324146 

- 0.3367190 

0.0342828 

1.0000 

- 1.0000000 

- 0.0063047 

- 1.0000000 

0.0353790 

- 1.0000000 

0.0351764 


Table 2: Numerical solutions obtained along the horizontal centreline. 



Re = 

100 

Re = 

1000 

Re = 

10000 

X 

V 

p 

V 

p 

V 

p 

0.0000 

0.0000000 

0.0194764 

0.0000000 

0.0422581 

0.0000000 

0.0458727 

0.0312 

- 0.0522265 

0.0225569 

- 0.1460808 

0.0432335 

- 0.3915446 

0.0441632 

0.0391 

- 0.0648647 

0.0232348 

- 0.1877457 

0.0432677 

- 0.4113208 

0.0425346 

0.0469 

- 0.0770068 

0.0238568 

- 0.2277794 

0.0431372 

- 0.3964249 

0.0408660 

0.0547 

- 0.0887644 

0.0244285 

- 0.2651327 

0.0428066 

- 0.3704853 

0.0393322 

0.0937 

- 0.1401717 

0.0264122 

- 0.3742712 

0.0378014 

- 0.3110442 

0.0330918 

0.1406 

- 0.1809713 

0.0265146 

- 0.3396911 

0.0280398 

- 0.2709560 

0.0261939 

0.1953 

- 0.1951736 

0.0233699 

- 0.2491894 

0.0186821 

- 0.2239381 

0.0189608 

0.5000 

0.0501305 

0.0000000 

0.0276485 

0.0000000 

0.0077100 

0.0000000 

0.7656 

0.1382907 

0.0082429 

0.2488983 

0.0295021 

0.2001465 

0.0187390 

0.7734 

0.1375760 

0.0084751 

0.2527025 

0.0309354 

0.2059849 

0.0197800 

0.8437 

0.1216119 

0.0102571 

0.2558328 

0.0427657 

0.2596504 

0.0302368 

0.9062 

0.0900690 

0.0115459 

0.2121234 

0.0493963 

0.3135677 

0.0408176 

0.9219 

0.0789643 

0.0118740 

0.1947518 

0.0503543 

0.3254261 

0.0435931 

0.9297 

0.0729086 

0.0120435 

0.1846614 

0.0507312 

0.3275184 

0.0449388 

0.9375 

0.0664792 

0.0122188 

0.1732839 

0.0510457 

0.3253029 

0.0462217 

1.0000 

0.0000000 

0.0139673 

0.0000000 

0.0518188 

0.0000000 

0.0508029 
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Figure 7: Particle trajectories: Re=100 (left), Re=1000 (middle) and Re=10000 (right). 


Figure 7 visualises the flow field by means of particle trajectories. The vortex dynamics is similar to 
the classic lid-driven cavity case, except that there is a small tertiary vortex at the top left corner for 
Re=10000, in agreement with [45]. 

Figure 8 shows the truncation errors (calculated via (18)) of each equation on the 1024x1024 grid 
for the various Reynolds numbers. These images will be helpful later on in order to interpret why each 
refinement criterion chooses to refine the specific areas that it does. 

6.2. Assessment of the efficiency of local refinement in reducing the discretisation error 

To assess the various local refinement strategies, some measure of the discretisation error must be 
plotted against the number of CVs of the grid in each case. The best strategy would be that which achieves 
the greatest discretisation error reduction with the smallest number of CVs. When the discretisation error 
is plotted against the number of CVs, the curve with the steepest slope corresponds to the most efficient 
refinement strategy. Three kinds of such plots are included in Section 7: 

1. The discretisation error is plotted against the total number of CVs at a point where the grid is 
locally dense. 

2. The same as above, but for a point where the grid is locally coarse. 

3. The norm || || i of the whole discretisation error in the domain is plotted against the number of 
CVs. This norm is defined as: 


|| l 


1 K 

1=1 


(32) 


where D is the total volume of the computational domain, and and are the discretisation error 

and volume of CV i of grid h. 

In case (1) the discretisation error reduction is partly or mostly due to grid refinement in the vicinity 
of the particular point. In case (2) the reduction is mostly due to grid refinement at distant regions from 
where the discretisation error is transported to that particular point. In case (3) the effect on the whole 
domain is investigated. Depending on the actual application, one may either be interested in achieving 
high accuracy in a limited region of interest, or in globally achieving high accuracy throughout the domain. 

For cases (1) and (2) the points are selected among those listed in Tables 1 and 2, and the discretisation 
error is calculated based on the values displayed in the tables. For case (3) things are more complex because 
the “exact” solution must be calculated at the centre of each CV of the locally refined grid, as required 
by equation (32). The following approach was used: First, the 2048x2048 solution was restricted to the 
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Figure 8: Base-10 logarithm of the absolute value of the truncation errors of the x —momentum equation (left), of the 
y —momentum equation (middle) and of the continuity equation (right), for Re=100 (top), Re=1000 (middle) and Re= 10000 
(bottom). 
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Figure 9: Interpolation stencils. Possible interpolation locations are marked with x. CV centres or boundary face centres 
or corner points of the 1024x1024 grid are marked with •. The CV which contains the interpolation point is shown shaded. 
See text for more details. 


1024x1024 grid using a third-order accurate restriction operator which will be described below. This 
restricted solution together with the 1024x1024 solution were used to obtain the “exact” solution on grid 
1024 xl024 via Richardson extrapolation assuming second-order convergence. The solution obtained on 
each composite grid was compared against this “exact” solution to calculate the discretisation error on 
each CV of the composite grid. Since on cell-centred grids, such as those used in the present study, the CV 
centres of different levels do not coincide, the “exact” values at the centres of the CVs of the composite 
grids had to be interpolated from the values at the 1024x1024 grid in general. The error introduced by 
this interpolation should be smaller than the discretisation error, otherwise the interpolated values cannot 
be considered as “exact”. Therefore, since a second-order accurate finite volume method was used here, 
a third-order accurate interpolation operator was chosen which will now be described. 

Suppose the centre of a CV of a composite grid, where the “exact” solution is to be interpolated 
(one of the x points shown in Figure 9). First, a CV of the 1024x1024 grid which contains the given 
point is identified (shown shaded in Figure 9). Actually, the following three situations are possible: If 
the composite grid CV is of a level coarser than the finest level of the 1024x1024 grid then its centre will 
coincide with one of the corners of the containing CV (such a point also has alternative container CVs); 
if it is of the same level as the finest level of the 1024x1024 grid then its centre will coincide with the 
centre of the containing CV; and if it is of a finer level (equal to the finest level of the 2048x2048 grid) 
then it will lie in the interior of the containing CV. The number of CVs of the composite grid may be 
large, but locating the container CV for each of them is not expensive because it is performed using a tree 
algorithm: First a container CV of the coarsest level is found, then the container CV of the immediately 
finer level is sought among that CVs children, and so on until a container CV of the 1024x1024 level is 
obtained. 

The neighbours of the containing CV are also identified and their centres are recorded (shown as black 
dots in the left part of Figure 9). If the containing CV is a boundary CV or a corner CV then some of 
the neighbours will be missing, in which case the corresponding boundary face centres or corner point 
are selected instead (also marked as black dots in Figure 9 middle and right parts respectively). In any 
case, we have recorded the centre of the containing CV (£ 1 , 2 / 1 ) and the 8 neighbouring points (^ 2 ,^/ 2 ) 
to (£ 9 , 2 / 9 ). The “exact” values 0 i, 02 , • • • ,09 are also known at these points, where 0 is any of the flow 
variables - either by Richardson extrapolation or from the boundary conditions. Suppose also that (£ 0 , Vo) 
is the point where we want to interpolate the function. 

The exact solution is reconstructed from these values in the vicinity of these points by a quadratic 
polynomial of the form: 

4>(x,y) = ci + c 2 x + c 3 y + c 4 x 2 + c 5 xy + c 6 y 2 (33) 

where x = x — xo and y = y — 2/0 • The coefficients q are selected such that the values of the polynomial 
at the 9 neighbouring points are as close as possible to the “exact “ values 0 1 , 02 etc. Therefore the 
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coefficients q should satisfy the following linear system A • c = b if possible: 
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This system has 9 equations but only 6 unknowns, and is therefore over deter mined. It is therefore 
not possible in general for the polynomial to return exact values at the 9 neighbouring points, but it is 
possible to minimise the error by solving the above system in the least squares sense, i.e. to find c for 
which the norm ||6 — Ac \\2 is minimised. This is achieved by computing the reduced QR factorisation of 
A, A = QR (we used the modified Gram-Schmidt algorithm) and then solving the 6x6 upper triangular 
system Rc = Q T b for c. For more details see [46]. The relative coordinates x, y are used rather than the 
actual coordinates x, y , because otherwise the column vectors of A are nearly parallel as the 8 points are 
very close to each other, and the QR factorisation is not accurate. This interpolation operator introduces 
an error of the order 0(/i 3 ) as we have verified with simple numerical experiments (not included). 

This interpolation operator was also used for the restriction of the 2048 x 2048 solution to the 1024 
x 1024 grid, in order to perform the Richardson extrapolation. In this case, the roles of the 2048 x 2048 
and 1024 x 1024 grids were interchanged, and the grid shown in Figure 9 would be the 2048 x 2048 grid. 
The centres of the 1024 x 1024 CVs would be those points marked with x that lie at the corners of the 
CVs of the 2048 x 2048 grid. Therefore, each 1024 x 1024 CV centre has 4 possible container CVs on 
the 2048 x 2048 grid. For maximum accuracy, the restriction of the variables was repeated for each of 
these 4 possible container CVs, and the “exact” value was taken as the average of the four interpolated 
values. 

6.3. Assessment of the effect of local refinement on the discretisation error distribution 

Even if local refinement is capable of efficiently reducing the mean discretisation error this does not 
necessarily mean that this error may not be large locally. Therefore, another useful property of local 
refinement would be its ability to smooth out the discretisation error distribution. Figure 8 shows that 
on structured grids there are huge variations in truncation error, by several orders of magnitude. Local 
refinement focuses on regions of high truncation error in an effort to reduce it and therefore the truncation 
error distributions become smoother on composite grids. It would be interesting to examine if a similar 
effect takes place as far as the discretisation error is concerned. 

Discretisation error distributions cannot be compared directly but they have to be normalised first by 
dividing them by their mean (given by (32)). So we define the normalised discretisation error as: 

e h = Tj —tt ' e h (35) 

\\C-h || 1 

The main property of which is of interest here is that the integral of its absolute value in the domain 
is independent of the order of magnitude of the actual discretisation error: 
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due to (32). In other words if we draw the e^(x,y) distribution as a surface in 3D space, the volume 
enclosed between this surface and the x — y plane always equals the area D of the domain (which equals 1 
in the present test case). This makes distributions comparable to each other not only among different 
refinement schemes but also between different Reynolds numbers and different flow variables (iq v, p ). 

Different e* h distributions are compared graphically by plotting dA ^ ^ against e* (see Figures 21, 22 
and 23, which are discussed in the next Section). Here A(e*) is the total area of the domain where the 
absolute value of the normalised discretisation error is less than or equal to e*. To calculate dA ^ - , first 
the interval [ 0 , e^ a J is split into a number of bins of size 0.1 each, e^ ax being the maximum absolute 
value of normalised discretisation error in the domain (bin 1 is the interval [ 0 , 0 . 1 ), bin 2 is the interval 
[ 0 . 1 , 0 . 2 ) etc.). Then the area of each CV of the grid is added to the appropriate bin according to the 

is calculated at the centre 


dA(e *) 
de* 


absolute value of the normalised discretisation error at that CV. Then 
of each bin as the total area of the bin divided by 0.1 (the width of the bin). 

The area of the domain where the absolute value of the normalised discretisation error is in the range 
[ 61 , 62 ] equals J^ 2 dA ^ ^ de*. In a graph such as that of Figure 21 that would equal the area below the 
curve and above the 6 -axis, between 61 and 62 . The area below the whole curve equals the area of the 
whole domain (/q 77 ™* ) de* = Q = 1. For a different test case with D % 1 one can use ^ dA } e * - instead 

of^). 

The mean value of |e*| in the domain, ^ f ^ |e*| dfl, equals 1 as can be seen from (36). If 6 * is relatively 
uniform within the domain then its distribution would be concentrated around e* = 1 . Otherwise, if it 
exhibits large variations, then its distribution would be more spread out. 

Also, the variation of |e*| in each case is quantitatively assessed using certain metrics (Tables 4, 5 
and 6 ). The first metric is e^ ax . However, this maximum value may be just an outlier and not be 
representative of the variation of the distribution. Therefore the 99-th percentile e* g% is also used. This 
is defined here as the upper endpoint of the bin i which is the lowest bin such that the sum of the areas of 
all bins from 1 to i is greater than or equal to 99% of the total area of the domain. In other words, e* g% 
is roughly such that the error |e*| is less than e gg% in 99% of the domain. Finally, the standard deviation 
is also used to quantify the variation of 6 *: 


a = 


\ 




2=1 


N 




(37) 


2=1 


where e^ = 1 is the mean value of |e*| in the domain. 


7. Results and discussion 


The results of testing eight different local refinement schemes for each Reynolds number are presented 
in this Section. These eight schemes will be referred to under the following descriptions: 

1 Q1 20°/ 0 XY c 

2 Q1 30°/ 0 XYC c 

3 Q2 20°/ 0 XY c 

4 Q2 20°/ 0 XYC c 

5 Q2 30°/ 0 XYC c 

6 Q3 20°/ 0 XY c 

7 Q2 20°/ 0 XYC n 

8 Q2 20°/ 0 XYC a 

The descriptions are interpreted as follows: Ql, Q 2 and Q3 denote the three local refinement criteria 
defined in Section 5.2 (Ql: truncation error; Q 2 : truncation error times CV volume; Q3: truncation error 
divided by main diagonal coefficient). Next comes the percentage R% of CVs that will be refined. That 
is, as described in Section 5.1, when a grid is to be refined the criterion quantity is calculated for every 
CV, the CVs are sorted in descending order according to the magnitude of the criterion quantity, and the 
first R% of this list are refined. This is repeated for all selected equations, which are stated next in the 
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above descriptions: XY means that only the momentum equations are used, while XYC means that the 
continuity equation is also used. The last character given in the descriptions describes the treatment at 
grid level interfaces: “c” means that only coarse CVs may be refined, “n” means that no refinement is 
allowed at all up to 4 CVs away from such interfaces, and “a” means that all CVs may be refined, i.e. 
there are no restrictions. See Section 5.3 for more details. 

In every case, the problem was first solved on a 32 x 32 grid (shown in Figure 6). Then the truncation 
error was calculated, and the grid was refined according to the selected refinement scheme to produce 
a composite grid. Then the refinement cycle was repeated on this composite grid (solution, truncation 
error estimate, refinement) to produce yet another composite grid. In fact this procedure was repeated 
six times. In every case the CVs of the finest level of the last (sixth) composite grid had the same size as 
those of the 2048 x 2048 structured grid. 

Six successive composite grids, times eight refinement schemes, times three Reynolds numbers produces 
a large total number of composite grids which are impossible to display here. Some of the grids are shown 
in Figures 10, 11 and 12. In every case the finest of the series of six grids is shown, and in fact it is the 
underlying grid 2 h which is shown, for clarity. 

Figures 13 - 15 show the convergence of ||e u ||i, the norm (32) of the discretisation error of u , with 
grid refinement for the selected Reynolds numbers. The graphs for the other variables (v, p) are similar, 
but in general mostly results on u are presented here, due to space limitations. Results for the other 
variables are included when appropriate. Figure 16 shows similar graphs but for the absolute value of the 
discretisation error of u at specific points. For each Reynolds number two points were selected among 
those listed in Table 1: One where the grids are coarse in general (left graph) and one where the grids are 
fine in general (the top-most point in fact - right graph). The results of Figure 16 must be interpreted 
with some caution because the reduction of the discretisation error at a specific point during a refinement 
cycle may be strongly influenced by whether the grid was refined or not at that particular location during 
that cycle. 

Table 3 contains results on the norms of the discretisation errors of the other variables as well, for a 
limited number of cases. These cases differ only in the quantity used as a refinement criterion (Ql, Q2 or 
Q3) while the other parameters of the refinement schemes are kept the same, to allow a direct comparison 
among these quantities. 

To investigate the effect of local refinement on the distribution of discretisation error across the domain, 
Figures 17, 18 and 19 depict the normalised ^-errors e* across the domain for Re = 100, 1000 and 10000 
respectively, for various cases. Also, Figure 20 shows the normalised discretisation error of pressure e* 
on the structured grids (top row) and on composite grids Q2 20% XYC c (bottom row) for the range of 
Reynolds numbers. Then, Figures 21 - 23 show the statistical distributions of the ^-errors, and Tables 4 
- 6 contain associated metrics. 

From this large amount of data, the following observations can be made concerning the influence of 
various parameters on the efficiency of local refinement: 

Effect of Reynolds number 

Figure 8 shows that as the Reynolds number increases, the variation of truncation error becomes 
greater: Both the regions of high r and the regions of low r become larger. The truncation error varies 


Table 3: Error norms on the finest grid for various cases. 
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(a) Ql, Re = 100 


(b) Q2, Re = 100 


(c) Q3, Re = 100 



1000 


(d) Ql, Re 



(g) Ql, Re = 10000 



(e) Q2, Re = 1000 



(h) Q2, Re = 10000 



(f) Q3, Re = 1000 



(i) Q3, Re = 10000 


Figure 10: The finest grids for cases: Ql 20% XYc (left column), Q2 20% XYc (middle column), and Q3 20% XYc (right 
column) (actually the underlying grids 2 h are shown). Re = 100 (top), Re = 1000 (middle), Re = 10000 (bottom). Different 
colours help to distinguish between levels. 
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(b) Re 


(c) Re = 10000 


(a) Re = 100 


Figure 11: The finest grids for case Q2 20% XYCc (actually the underlying grids 2 h are shown). Re =100 (left), Re = 1000 
(middle), Re = 10000 (right). 




(a) Q2 20% XYC n, Re = 1000 (b) Q2 20% XYC a, Re = 1000 

Figure 12: The finest grids for cases Q2 20% XYCn (left) and Q2 20% XYCa (right), Re = 1000 (actually the underlying 
grids are shown). 

by several orders of magnitude. For Re = 100 it is high near the lid and low everywhere else. As the 
Reynolds number increases, it becomes high also at the outer part of the main vortex, while in the interior 
of the vortex it becomes very small. These variations of r reflect on the variations of grid density as seen 
in Figures 10 and 11. 

Figure 13 shows that for Re=100 the benefits from local refinement are minimal concerning the total 
discretisation error in the domain, no matter what refinement scheme is used (Q 2 30% XYC c appears to 
be the most efficient - some schemes even appear to have a negative impact on the efficiency). However, 
Figure 16 shows that local refinement can be quite efficient at reducing the discretisation error locally 
where the grid is dense. The efficiency is sometimes very poor at coarse regions of the grid as also shown 
in Figure 16. Figure 17 shows that on structured grids there are very high discretisation errors near the 
lid, which are greatly reduced by local refinement. These errors concern only the u variable and are quite 
high relative to the mean error (Table 4: (e*) 99 % = 5.8 on the structured grid, (e*) 99 % ~ 3.5 on composite 
grids). On the other hand away from the lid the errors are not smoothed out, but, as shown in Figure 17, 
there are regions where e* is high no matter what refinement scheme is used. This is an indication of the 
complexity of the dependence of e on r (equation ( 2 )), which cannot be tackled by the simple refinement 
criteria tested here. Figure 21 shows that indeed there is some smoothing of e* caused by local refinement, 
as dA/de * is higher for the composite grids than for the structured grid near e* = 1 , and smaller for the 
composite grids than for the structured grid near e* = 0 and for large values of e*. But the difference is 


24 
















































































































































































































































































































































































































































Re = 100: 


\Eu\ 1 



Number of Control Volumes 


Figure 13: Convergence of the ||.||i norm (32) of the u —discretisation error with grid refinement for Re = 100. 


Re = 1000: \Eu\ 1 



Number of Control Volumes 


Figure 14: As for Figure 13 but for Re = 1000. 
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Re = 10000: 


|Eu| x 



Figure 15: As for Figure 13 but for Re = 10000. 


Table 4: Statistical indices of the distributions of normalised error, for Re=100. 


Re = 100 

(ej)max 

( e S)99% 

< 

(^)max 

( e S)99% 

< 

( e p)max 

( e J)99% 

(7 p 

Structured 

15.8 

5.8 

1.293 

5.52 

3.4 

0.897 

244 

3.1 

3.851 

Q1 20% XY c 

4.27 

3.6 

0.932 

3.85 

3.7 

1.015 

12.5 

4.0 

0.793 

Q1 30% XYCc 

4.00 

3.5 

0.878 

4.90 

3.8 

0.972 

34.3 

7.9 

1.622 

Q2 20% XY c 

4.73 

3.5 

0.873 

4.76 

3.5 

0.895 

13.0 

3.3 

0.585 

Q2 20% XYCc 

4.31 

4.0 

0.901 

4.30 

3.4 

0.875 

49.3 

6.1 

1.465 

Q2 30% XYCc 

5.01 

3.7 

0.868 

5.46 

3.7 

0.854 

66.2 

8.9 

2.150 

Q3 20% XY c 

4.76 

3.5 

0.870 

4.83 

3.5 

0.895 

33.1 

6.6 

1.331 

Q2 20% XYCn 

5.23 

3.9 

0.933 

5.39 

4.2 

0.966 

7.75 

3.8 

0.786 

Q2 20% XYCa 

6.41 

3.6 

0.874 

7.62 

3.9 

0.876 

51.3 

9.1 

1.787 


not dramatic. This mild smoothing is also reflected in the value of cr* in Table 4 1.3 on structured 

grid, ~ 0.87 on composite grids). No smoothing takes place for v though, as seen in the same Table. This 
is because there are no high e v regions near the lid. 

The pressure variable also deserves some attention for Re=100. The corresponding results of Table 
4 are difficult to interpret, but it must be taken into account that due to linear extrapolation to the 
boundaries the interpolation operator (33) may loose its accuracy and large errors may appear there. 
This is indeed the case since Table 4 lists a value of (e*) max = 244 but merely a value of (e*)gg% = 3.1 
which excludes the high errors at the boundaries. On the other hand, concerning the composite grids, 
relatively high errors occur also at level interfaces as may be seen in Figure 20, something which is not 
observed with u or v. In fact it is also not observed with p at higher Reynolds numbers. It is not known 
whether these high p errors have a significant effect on the overall efficiency of local refinement. 

As the Reynolds number increases, the efficiency of local refinement also increases - Figures 14 - 15. 
As mentioned in Section 2, the discretisation error is transported by convection and diffusion with the 
truncation error acting as source. Therefore as the Reynolds number increases and convection becomes 
stronger, the discretisation error is transported further away from its source regions of high truncation 
error, because convection tends to transport a quantity without reducing its magnitude, while diffusion 
attenuates the quantity that it transports. This means that reducing the source (truncation error) by 
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Re=100: Absolute U-error at x=0.5, y=0.1016 


Re=100: Absolute U-error at x=0.5, y=0.9766 




Number of Control Volumes 


Number of Control Volumes 


Re=1000: Absolute U-error at x=0.5, y=0.2813 


Re=1000: Absolute U-error at x=0.5, y=0.9766 




Number of Control Volumes 


Number of Control Volumes 


Re=10000: Absolute U-error at x=0.5, y=0.2813 


Re=10000: Absolute U-error at x=0.5, y=0.9766 




Number of Control Volumes 


Number of Control Volumes 


Figure 16: Convergence of |e u | at specific points, for various cases. The left and right columns show results at points generally 
located at coarse and fine regions of the composite grids respectively. 
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(d) Q1 20% XY c (e) Q2 20% XY c (f) Q3 20% XY c 


Figure 17: e* for various cases, Re =s 100. 


Table 5: Statistical indices of the distributions of normalised error, for Re=1000. 


Re = 1000 

(ej)max 

( € u) 99% 

< 

(^)max 

( e S)99% 

< 

( e p)max 

( e p)99% 

(7 p 

Structured 

10.5 

4.8 

1.007 

4.56 

3.1 

0.811 

10.6 

2.7 

0.568 

Ql 20% XY c 

3.61 

3.2 

0.814 

3.28 

3.0 

0.762 

3.81 

2.7 

0.548 

Ql 30% XYCc 

3.88 

3.6 

0.840 

4.03 

3.1 

0.783 

4.61 

2.7 

0.567 

Q2 20% XY c 

4.42 

3.6 

0.851 

3.83 

3.4 

0.816 

5.03 

2.8 

0.639 

Q2 20% XYCc 

4.57 

3.5 

0.855 

4.18 

3.5 

0.836 

5.67 

3.2 

0.649 

Q2 30% XYCc 

4.48 

3.7 

0.880 

3.84 

3.3 

0.808 

9.61 

2.7 

0.657 

Q3 20% XY c 

4.44 

3.6 

0.853 

3.83 

3.3 

0.819 

4.88 

3.0 

0.633 

Q2 20% XYCn 

7.02 

4.7 

1.053 

5.08 

4.0 

0.861 

3.49 

3.0 

0.644 

Q2 20% XYCa 

4.44 

3.4 

0.866 

4.29 

3.2 

0.839 

8.58 

3.0 

0.644 


local refinement will have a larger impact on the discretisation error of the whole domain when the 
Reynolds number is larger. This is also seen in Figure 16, which shows that for Re=1000 and especially 
for Re= 10000 the discretisation error reduces faster with local refinement even at coarse regions of the 
grid. Figures 18, 19 and 20 show that the numerical solution underestimates the strength of the main 
vortex on all grids. Local refinement again makes e* a bit smoother due its reduction near the lid, but e* 
and e* are not made smoother, e* in particular is a lot smoother even on structured grids as the Reynolds 
number increases. 

Effect of refinement quantity 

To compare the criteria Ql, Q2 and Q3 it is best to compare the results for schemes Q1 20% XY 
c, Q2 20% XY c and Q3 20% XY c. A decrease in CV volume reduces Q2 and Q3 both directly and 
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(d) Q1 20% XY c (e) Q2 20% XY c (f) Q3 20% XY c 


Figure 18: e* for various cases, Re = 1000. 


indirectly (through the reduction of the truncation error), but it only reduces Q1 indirectly. The result, 
as Figure 10 shows, is that Q1 produces grids that are very fine at regions oh high r and very coarse in 
regions of low r, greatly reducing e* near the lid but increasing it further away. Figure 19 shows that Q1 
20% XY c is the only scheme which allows the development of errors in the middle of the domain, due to 
the coarseness of the grid there. Q2 and Q3 produce smoother grids and Figures 13 - 15 show that this 
makes them more efficient than Ql. An exception is Q1 20% XY c for Re = 10000 which although it is 
less efficient than Q2 and Q3 at coarser grids, it catches up in efficiency as the grids become finer, and 
in fact it has the steepest slope suggesting that if even more refinement cycles were carried out it could 
surpass the other quantities in efficiency. Figure 16 shows that Ql may be quite efficient at particular 
fine-grid points. 

Figure 10 shows that Q2 and Q3 produce very similar grids, as expected according to the discussion 
in Section 5.2, since all CVs of the grid have the same aspect ratio. The differences in the grids they 
produce are minimal at Re = 100 but are slightly more noticable at Re = 10000, which is not surprising 
since the convection terms of the main diagonal coefficients are stronger in the latter case. Figures 15 
and 16 show that for Re = 10000 Q2 is slightly more efficient which is interesting, since it suggests that 
the activation of the convection terms of Ap^p in quantity 3, which reduces the chance of a CV being 
refined if the velocity is high, makes things worse rather than better. Figures 13, 14 and 16 show that 
both quantities produce very similar results for Re = 100 and 1000. 

Effect of continuity equation 

To investigate the effect of whether the continuity equation is used in the refinement scheme we 
compare the schemes Q2 20% XY c and Q2 20% XYC c. Figures 13 - 15 reveal that the use of the 
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(d) Q1 20% XY c (e) Q2 20% XY c (f) Q3 20% XY c 


Figure 19: e,* for various cases, Re = 10000. 


Table 6: Statistical indices of the distributions of normalised error, for Re=10000. 


Re = 10000 

(ej)max 

( e £) 99% 

a u 

(^)max 

( € v)99% 

(T v 

( e p)max 

( e S)99% 

fj p 

Structured 

9.80 

3.5 

0.901 

6.97 

3.5 

0.862 

4.01 

2.6 

0.592 

Q1 20% XY c 

5.62 

3.8 

0.875 

6.24 

3.9 

0.858 

1.67 

1.6 

0.395 

Q1 30% XYCc 

3.19 

2.5 

0.758 

3.15 

2.5 

0.753 

2.03 

1.8 

0.385 

Q2 20% XY c 

8.91 

3.2 

0.779 

3.56 

3.1 

0.757 

2.97 

2.2 

0.531 

Q2 20% XYCc 

3.80 

2.8 

0.758 

3.07 

2.8 

0.733 

2.93 

2.3 

0.531 

Q2 30% XYCc 

5.83 

2.9 

0.778 

4.21 

2.7 

0.751 

3.00 

2.2 

0.521 

Q3 20% XY c 

4.34 

3.1 

0.754 

3.40 

2.8 

0.720 

2.67 

2.1 

0.474 

Q2 20% XYCn 

3.72 

3.0 

0.817 

3.68 

2.8 

0.767 

2.16 

2.1 

0.483 

Q2 20% XYCa 

5.38 

3.2 

0.804 

4.25 

3.2 

0.803 

3.38 

2.4 

0.555 
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(d) Re=100: Q2 20% XYCc (e) Re=1000: Q2 20% XYCc (f) Re=10000: Q2 20% XYCc 


Figure 20: e* for various cases. 


continuity equation makes local refinement slightly more efficient for Re = 100 and slightly less efficient 
for Re = 10000, compared to the case where the refinement decision is based only on the momentum 
equations. For Re = 1000 the difference in efficiency is marginal. The increase in the number of CVs of 
the grid caused by the inclusion of the continuity equation in the refinement criterion is 36% for Re = 
100, 28% for Re = 1000 and 31% for Re = 10000 (compare also Figure 11 against the middle column of 
Figure 10). Also noteworthy is that Table 4 suggests that the inclusion of the continuity equation worsens 
the problem of high e* at level interfaces. 

Effect of refinement ratio 

Comparing schemes Q2 20% XYC c and Q2 30% XYC c in Figure 13 reveals that increasing the ratio 
form 20% to 30% improves the efficiency slightly, just like including the continuity equation which also 
increases the number of CVs of the grid. A much greater improvement is observed among schemes Q1 
20% XY c and Q1 30% XYC c which includes the combined effect of including the continuity equation 
and increasing the refinement ratio. For Re = 1000 there appears to be no benefit from increasing the 
refinement ratio as far as quantity Q2 is concerned, but there is an important benefit combined with the 
use of the continuity equation for Ql. For Re = 10000 the increase of the refinement ratio reduces the 
efficiency, just like including the continuity equation. This shows that for high Reynolds numbers it is 
more efficient to restrict local refinement to small regions of high momentum truncation error. 

Effect of level interface treatment 

Comparing schemes Q2 20% XYCc, Q2 20% XYCn and Q2 20% XYCa it is clear that it is most 
efficient to allow refinement only at the coarse side of a level interface. Not allowing refinement at ah 
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Figure 21: Distributions of absolute normalised u-error on the finest grids, Re = 100. 
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Figure 22: Distributions of absolute normalised ^-error on the finest grids, Re = 1000. 
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Figure 23: Distributions of absolute normalised u-error on the finest grids, Re = 10000. 


produces terrible results at all Reynolds numbers, as Figures 13 - 15 show. The reason is clear from 
Figure 12(a): Refinement is prevented from occurring in large areas of the domain. On the other hand, 
not posing any restrictions at all produces grids like the one shown in Figure 12(b) where level interfaces 
themselves constitute refinement nuclei and consume many of the available CV refinements. This is 
expected to become worse as refinement cycles progress, because the truncation error reduces at the rest 
of the domain due to refinement, but not at level interfaces where r E 0(1). An advantage of criteria Q2 
and Q3 in this respect is that the size of the CV alone plays a role in the selection of the CVs to be refined. 
Therefore even though the truncation error does not decrease with grid refinement at level interfaces, the 
quantities Q2 and Q3 do decrease, thereby decreasing the likelihood of refinement in these interface CVs. 
Figures 13 - 15 show the rather surprising result that scheme Q2 20% XYCa is only slightly less efficient 
than Q2 20% XYCc. The e* distributions in space of Q2 20% XYCa are very similar to those of the 
other refinement schemes as shown in Figures 17 - 19 only that they are rugged, but the ruggedness is 
smoothed out as the Reynolds number increases. 

8. Conclusions 

For the present study a number of numerical experiments were conducted involving the solution of the 
regularised lid-driven cavity problem using various local grid refinement schemes for a range of Reynolds 
numbers. All local refinement schemes are based on an estimate of the truncation error.The goal was 
to assess the effect of varying various local refinement parameters on the efficiency of the finite volume 
method. To this end the regularised lid-driven cavity problems were solved to a great level of accuracy. 
Several conclusions can be drawn from the results presented in this paper. 

First of all, truncation error-based local refinement appears to offer little benefit on a global scale when 
the Reynolds number is low, although it can be efficient at a local scale. However, local refinement becomes 
quite efficient as the Reynolds number increases. It is suggested that this is due to the discretisation error 
being convected farther away from the source locations (of high truncation error) at higher Reynolds 
numbers, so reducing the sources causes a greater reduction in the global discretisation error. 

Refinement criteria that are proportional to the CV volume in addition to the truncation error are 
found to be more efficient than simply using the truncation error by itself. Using the product of the 
truncation error times the CV volume has been found to be better than using the truncation error 
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divided by the main diagonal coefficient of the linearised operator. 

Using the truncation error of the continuity equation in addition to those of the momentum equations, 
and increasing the percentage of CVs that are refined in each refinement cycle increases the efficiency 
when the Reynolds number is low, but decreases it when it is large. 

Despite the fact that the present local refinement schemes attempt to smooth out the truncation error 
distribution in the domain, it was observed that they only slightly smooth the distribution of discretisation 
error. This is an indication of the complexity of the r — e relationship. 

Also, an effort was made to deal with the issue that at level interfaces r G 0(1). Although the exact 
effect of these high truncation errors on the discretisation error is not completely clear, the present results 
show that this issue is best dealt with by allowing refinement to occur only at the coarse side of a level 
interface rather than either not allowing refinement at all there or not imposing any restrictions. However, 
in the latter case of not restricting refinement at level interfaces, rather surprisingly, the efficiency does 
not appear to suffer severely by the pile-up of refinement levels on top of the original level interface. 

Finally, we would like to make some general comments on the refinement strategy adopted and to 
discuss some implications in the case that a transient problem was solved instead. Most of the previous 
researchers that based their refinement criterion on either the truncation error or the residual used a 
threshold to select the CVs to be refined; a CV is refined if the refinement quantity there is above the 
threshold. The refinement cycles continue until the refinement quantity is below the threshold at all CVs, 
or until a maximum allowed number of levels has been reached. We did not follow this procedure because 
of two reasons: The first is that the truncation error does not decrease at certain regions of the grid 
such as at level interfaces and at domain boundaries, and therefore it is impossible for the refinement 
criterion quantity to drop below the threshold everywhere. The second is that due to the complexity of 
r — e relationship it is very difficult to relate a desired accuracy to a certain threshold. Therefore, in 
our finite volume scheme such a procedure would cause the following: The refinement procedure would 
stop in the interior of grid levels when the quantity dropped below the threshold, but the refinement 
cycles would continue since the quantity would be above the threshold at level interfaces and at domain 
boundaries. Eventually the cycles would stop when the maximum allowed level was reached. One would 
then have to decide on choosing a particular maximum allowed level. By allowing more and more levels 
the number of CVs of the grid would continually increase (at level interfaces and domain boundaries), 
but the discretisation error would not converge to zero since CVs would not be refined in the interior of 
the levels where the refinement quantity is below the threshold. The discretisation error would converge 
to a certain non-zero value, as observed e.g. in [11]. 

On the other hand, the aim of the present study is to compare a number of refinement criterion 
quantities, and therefore the most natural approach would be to refine the CVs where this quantity is 
largest. The comparison would be most accurate if it happened on a CV-by-CV basis, i.e. if during a 
refinement cycle only one CV was refined, the one where the refinement quantity is largest. Then by 
observing how fast the discretisation error reduces as the number of refinement cycles increases, one could 
evaluate the various refinement quantities. Of course, such a procedure is completely impractical due to 
the overhead involved in each refinement cycle (the data structures have to change, the equations have to 
be solved on the new composite grid, the truncation error has to be calculated afresh etc.). However, the 
next best thing that one can do which is also practical, is to not refine a single CV in a refinement cycle 
but to refine many CVs selected on the basis that their refinement quantity is larger than at the rest of 
the CVs of the grid. So in the present study the number of these CVs was selected to be a percentage R 
of the total number of CVs of the grid. This allows a clearer comparison between the various refinement 
quantities than using a threshold, because each refinement quantity would have its own threshold and it 
is not straightforward to make these thresholds equivalent. Also, by allowing more and more refinement 
cycles to occur this procedure (which uses the percentage R ) causes the discretisation error to converge 
to zero, if the refinement quantity is proportional to the CV size, (as the present results verify), since 
there is no threshold below which a CV will not be refined. Furthermore, the percentage R combined 
with a chosen fixed number of refinement cycles makes the number of CVs of the final composite grid to 
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be approximately known a priori; it approximately equals K • (1 + 3 R) n where K is the number of CVs of 
the initial grid, and n the number of allowed refinement cycles (when a CV is refined it is replaced by 4 
child CVs, so the total number of CVs increases by 3), so we have two parameters, R and n which can be 
varied simultaneously to get the desired grid size. Therefore, the refinement strategy adopted here can be 
described as aiming to produce the best possible accuracy for a given selected grid size. This is achieved 
better if R is small, but n must not become so large that the overhead of the refinement procedure causes 
the overall efficiency to degrade. A disadvantage of the method is that it does not provide much control 
over the discretisation error, but the same also holds even when a threshold is used. 

The present study is concerned only with the steady-state case. Transient problems would have 
additional complexities. First of all the truncation error would contain also components which would 
be a function of the time step size At, in addition to those which are functions of the grid spacing h. 
The present truncation error estimation method which uses the coarse grid 2 h would only identify the 
part of the truncation error which is due solely to the spatial discretisation. In order that the whole 
truncation error be calculated the coarse grid operator should also use a time step of 2A t. In this case 
local refinement should also be performed in time, by using smaller time steps where appropriate. To 
reduce the complexity one can limit the local refinement to the spatial grid only, trying to reduce only 
the spatial components of the truncation error, but this procedure cannot be expected to reduce the 
discretisation error below a certain point, as the temporal components of the truncation error do not 
decrease. The refinement procedure described here could be used within each time step if an implicit 
temporal discretisation scheme is chosen, but not without modification. Indeed, using a refinement 
percentage R time step after time step would cause the grid to become huge in terms of number of CVs 
very soon. Since the aim of the present method is to achieve the best accuracy with a given grid size, the 
procedure could be adapted as follows: During the first time step an initial coarse grid would be adapted 
using selected values of R and n, so that the desired grid size is reached, and during the following time 
steps there would be as many CV refinements as there are unrefinements, so that the total number of CVs 
of the composite grid remains constant, with the aim that the refinement quantity is as low as possible in 
the domain. The increase in truncation error that an unrefinement would cause can be estimated from the 
truncation error estimate at the child CVs assuming r E 0(h p ). Therefore in this case it is necessary that 
the method has the ability to perform unrefinement, that is to remove 4 child CVs to recover their parent. 
This is an issue with its own complexities that is beyond the scope of the present work; implementation 
details can be found e.g. in [47, 48]. 
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